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1 INTRODUCTION 


Half (or more) of a flight vehicle’s drag arises from skin friction. It is therefore highly desirable to 
have accurate measurement tools that enable the measurement of skin friction distributions on surfaces. 
To date, various skin friction measurement methods have been developed. One is the fringe-imaging 
skin friction (FISF) technique, which was originally developed by Monson and Mateer (ref. 1) and has 
been used successfully in many different flows. With this technique, the flow pattern of an oil patch 
placed on the surface of the wind tunnel model is observed. To determine the skin friction magnitude, the 
thickness variation of the oil is measured using interferometry and related to the skin friction magnitude 
through lubrication theory. 

In Squire’s (ref. 2) theoretical work on the motion of a thin oil sheet under a steady boundary layer, 
a nonlinear partial differential equation that relates skin friction to the thickness variation of the oil was 
derived. Tanner and Blows’ (ref. 3) integral approach yields a linear equation for oil layer thickness as 
a function of skin friction, time, oil viscosity, and position and was shown to be the linearized solution 
to a reduced form of Squire’s equation. It forms the basis for the FISF technique. 

Zilliac (ref. 4) performed a numerical simulation of Squire’s reduced equation and proved the 
validity of the linearized solution. The numerical results are also in accordance with experimental data 
(ref. 4) although no comparisons are shown at an early time when surface tension effects are thought 
to be large. 

The main assumptions, which were made to obtain the linearized equation used in the FISF tech- 
nique, are that the oil flow is an incompressible, two-dimensional, slow viscous motion and that surface 
tension as well as wall adhesion effects can be neglected. The last two have been omitted entirely in 
previous studies. 

In this paper a modified version of the computer program RIPPLE (ref. 5) is used to address surface 
tension and wall adhesion effects. RIPPLE computes finite difference solutions for incompressible, two- 
dimensional, laminar viscous flows with free surfaces. The volume-of-fluid (VOF) method (ref. 6) is 
used for surface-tracking and the continuum surface force (CSF) model (ref. 7) for inclusion of the 
surface tension effects. RIPPLE’s capabilities have been shown in many examples; the collision of two 
water rods and the low-gravity jet-induced tank flow are two of them (ref. 5). 

In the current study, the influence of surface tension and wall adhesion on the accuracy of the FISF 
technique is investigated. In particular, the validity of the linear relation between constant shear stress, 
downstream position, oil viscosity, and the height of the oil film in the presence of surface tension is 
questioned. 



2 MODELING APPROACH 


In the following we will give an overview of the flow problem studied in this work. Previous 
approaches, which will be used for comparison with the obtained computational results, are discussed 
briefly. The physical model (in particular the phenomena of surface tension and wall adhesion) ^eluding 
the defining equations is illustrated, and the computational technique with special focus on the O 
and CSF model is presented. The latter is extended for tangential tensile forces arising from skin 

friction. 


2.1 Considerations 

Figure 1 shows a schematic representation of the flow being studied. A line of oil is placed on 
a flat plate. A nontransient, turbulent air boundary layer forces the oil to spread in the do— 
direction. The oil height is so small that it remains entirely in the viscous sublayer. The Realized 
drop shown in figure 1 is a hemispherical drop with the same fluid volume as the initial drop 
it is P solely used for the purpose of defining a suitable characteristic length. The initial drop is the oil 
distribution at time t = 0 for the RIPPLE calculations (and also for the Lax-Wendroff computations). 

The linear equation used in the FISF technique, results of Zilliac’s numerical a PP ro ^°^ m P^ g 
the thinning rate of the oil film, and limited experimental data will be compared with RIPPLE results. 
The linear solution and Zilliac’s Lax-Wendroff-based simulations will be discussed briefly. 

2.1.1 Linear Solution and Lax-Wendroff Simulation 

Squire (ref. 2) performed a creeping-flow-conservation-of-mass analysis of an oil film flow and 
obtained an expression for the thinning rate of a three-dimensional oil patch. For a two-dimensional 



Figure 1 . Oil drop underneath an air boundary layer. 
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case, this equation becomes 


duaijA h ? 1 h? dp 

dy ) y=h 2 pw 3 dx ^ 

where A is the ratio of the air viscosity and oil viscosity (A = Pw/ P q\])- Squire showed by order of 
magnitude considerations that the pressure gradient in the oil normal to the wall is negligible ( dp/dy = 
0), in other words: p(x ) a j r = p(x) 0 j]. Hence the pressure gradient dp/dx and the velocity gradient 
(diiw/ dy)y~f l in equation (1) can be obtained directly from boundary-layer theory. 

The boundary conditions at the air-oil interface (y = h) are: 

(1) that the air and oil velocities are the same and 

(2) that the shear stresses in the oil and the air are equal at the free surface. 


1 dh _ _d_ 
A dt dx 


Equation (1) can be further simplified by assuming that the pressure gradient term is small and by 
alignment of the coordinate system with the surface shear stress direction. This yields 


dh 1 _d_ (h? \ 

dt p 0 \\ dx \ 2 Twan J ^ 

with r wa u = {dua^l dy)y—} l . Equation (2) can be solved analytically resulting in the linear expression 


^ ^Moil 

^wall 


(3) 


In the FISF technique, h(t) is measured at a certain location x to obtain T wa u, using equation (3). This 
linear solution of the partial differential equation (2) assumes an infinite initial oil thickness (h = oo at 
t = 0) and a constant shear stress r wa u. (Note: in the above expression, x = 0 at the leading edge of 
the oil.) 


Zilliac (ref. 4) applied the Lax-Wendroff algorithm to solve a finite difference form of equation (2). 
The linear equation and the Lax-Wendroff results were nearly identical at later times (t > 4 s). However, 
at early times, a comparison cannot be made because of the infinite initial height condition required to 
obtain the linear solution. 


2.1.2 Dimensionless Numbers 

In flows where inertial, interfacial, gravitational, and viscous forces play a role, the following 
four dimensionless numbers are frequently defined: the Reynolds number. Re = uljv; the Bond 
number, B = pgfi/a; the capillary (or Taylor-Saffman) number, Ca = up/a; and the Weber number 
W* 2 = R- e • Ca, where l is the characteristic length-scale and a is the surface tension coefficient. 

The Bond number can be interpreted as a ratio of gravity and surface tension. Thus a high Bond 
number indicates that gravitational forces have a much higher impact on the flow than interfacial forces. 


3 



Similarly, the capillary number indicates the importance of viscous forces relative to surface tension, 
while the Weber number compares inertial forces and surface tension. 


Referring to figure 1, one can see that a characteristic length is hard to define because the geometry 
of the flow is not static (an initially defined characteristic length would change with time). Howe™ •» 
rough estimate of the four relevant dimensionless numbers can be made if the radius of an 'deal 
drop" is chosen as the characteristic length. Note that the velocity u in the above dimension ess 
numbers refers to the velocity of the advancing contact line, which is in the cunent study approximarely 
„ = 0 84 mrn/s, and not to a free stream velocity normally used for obtaining the Reynolds numbers. 
The oil viscosity is given with ife = 50 Cs. Given these assumptions, the Reynolds, Bond, capillary, 
and Weber numbers become: 

Re = = — = 0.0063 , 

viscous v 

B = gravitational = P9? = Q 063 t 
interfacial o 
_viscous_ = «/* = 0 0002 
interfacial a 


We = 


inertial 


= Re - Ca = 1.27 • 10' 


respectively. 

The values of Ca and We are well within the typical ranges for flows dominated by mterfacial and 
viscous forces (ref. 8). This implies that surface tension forces may have a significant influence on e 

oil flow locally. 


2.2 Defining Equations 

2.2.1 The Differential Navier-Stokes Equations in Conservation Form 

The flow field in the oil film is governed by the Navier-Stokes equations for two-dimensional, 
incompressible, and viscous flows, 

V • V = 0 ( - 4 - 


— + V- (VV) = --Vp + -V -r + g + ^-Fi ( 5 ) 

dt ' J P P P 

Here V = V(x,y,t ) is the velocity field, g is the gravity acceleration^ F b is a body force, p is the 
pressure, r is the viscous stress tensor for a Newtonian fluid, r = p[V • V + (V • V) 1 , and p is the oi 
density. 

2.2.2 Surface Tension 

Within a fluid, the intermolecular forces of attraction are balanced in all directions. However, for 
the fluid molecules at or near the surface, the cohesive forces are unbalanced due to the discontinui y 
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of fluid properties at the interface (fig. 2). The unbalanced forces of intermolecular attraction at the 
surface cause the interface to behave like it were an elastic membrane under uniform tension. The 
magnitude of the tensile force per unit length of a line on the interface is called surface tension a. The 
value of a depends on the fluids in contact and the temperature. In order to take the effects of surface 
tension into account, special boundary conditions must be applied to the Navier-Stokes equations at 
an interface. Assuming that the heat flux and mass transport through the surface are negligible, the 
boundary condition that must be satisfied at the boundary between the two fluids yields an equilibrium 
expression for all forces acting on both sides of a surface element at an instant in time. 


-l \ OCT 

^ (Pfluid CKj Tfluid * 71 “b — 71 .Pvapor ^vapor * 71 


( 6 ) 


fluid forces 


vapor forces 


where n is the free surface unit normal vector directed into the fluid, taum is the viscous stress tensor 
in the fluid, T vapor is the viscous stress tensor in the vapor, and n is the local curvature calculated from 
(ref. 9) 

« = -(V-n) (7) 

Notice that da/ ds in equation (6) is the gradient of the surface tension coefficient along the surface. 

Further, assuming that the stress tensor is continuous across the interface (r fiuid = r vapor at the 
surface) (ref. 10) and that a can be regarded as constant ( da/ds = 0), equation (6) reduces to 


P fluid Pvapor — Ap — OK (8) 

The pressure jump A p in the above expression has to be applied as a boundary condition of equation (5) 
at the fluid interface. 


2.2.3 Wall Adhesion and Contact Angle 

When a drop of liquid is placed on a solid surface, the liquid spreads out under the influence of 
gravity. The rate of spreading depends on the relative sizes of intermolecular cohesive and adhesive 
forces. For example, if the fluid molecules possess a greater affinity for each other than for the solid 
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(i.e., cohesive forces are greater than adhesive forces) the drop tends to deform slightly but it does not 
“wet” the solid. 

In the opposite case, where the adhesive force between the solid and the liquid is greater than 
the cohesive forces between the fluid molecules, the drop spreads out until it reaches a state of static 
equilibrium. At this point, the surface of the drop forms an angle of contact (6 S ) with the solid surface 
(fig 3) which can be measured experimentally (ref. 11). Young (ref. 11) assumed that the magnitude 
of the static contact angle depends on the surface tension coefficients of the three material boundaries 

and In order for the system to be in state equilibrium, the horizontal 

components must sum to zero at the contact line, 

O vapor/solid — ^liquid/solid — ^liquid/vapor cos @S = 0 ^ 

Although this equation ignores the force balance normal to the solid surface by assuming that the solid 
acts like a rigid body, it is a good approximation from the point of view of macroscopic capi ary 
concepts. In the case of a moving contact line, the resulting dynamic contact angle is different from 
the static contact angle; it is at the very least a function of 9 S and the capillary number. 

The general equations used to determine the static and dynamic contact angles are still a subject 
of research. Some success has been achieved using empirical methods but a general theory of dynamic 
contact lines has not been derived (ref. 11). 


2.3 Computational Technique 

RIPPLE is a computer program for simulating transient, two-dimensional, incompressible fluid 
flows governed by equations (4) and (5). Free surfaces are tracked using a VOF technique that was 
developed at the Los Alamos National Laboratory by Hirt and Nichols (ref. 6). Surface tension forces 
are incorporated applying the CSF model (ref. 9). 

23.1 VOF Method 

The VOF technique for representing free surfaces assigns a characteristic value F to each compu- 
tational cell that specifies whether the cell is a fluid cell (cells completely filled with fluid), * surface 
cell (cells partially filled with fluid), or a void cell (empty cells). In two-fluid problems, such as the 
present one, the gas phase is not fully modeled by RIPPLE, and it is referred to as the void region^ In 
RIPPLE, F = 1 for a fluid cell, F = 0 for a void cell, and 0 < F < 1 in cells containing the surface. 



Figure 3. Surface tension forces at the contact line. 
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In these cells, F represents the fractional volume of the cell occupied by fluid. The VOF function 
F = F(x,t) is a scalar field advected as a Lagrangian invariant in an Eulerian scheme. 


dF 

dt 


dF 8F dF n 

m + u te +v %; =0 


( 10 ) 


where F(x, 0) is known. Equation (10) has to be solved in addition to the flow governing equations (4) 
and (5). The solution algorithm for equation (10) has to assure that the discontinuous property of F 
(which is a fundamental premise in multi-fluid problems, where it is desirable to maintain a sharply 
defined interface) is preserved during the computation. Standard advection algorithms, which compute 
fluxes algebraically, tend to spread discontinuities owing to numerical diffusion. Even minimizing the 
diffusion with the objective of restricting the spreading to a width of two or three cells does not match 
the capabilities of VOF advection methods, which typically remain within one cell width. 


In the past two decades, several VOF algorithms have been developed to advect volume fractions. 
Theses methods all have in common that they compute the variation of F in a cell by reconstructing 
the surface from VOF data and computing the volume fluxes geometrically. The surface reconstruction 
method plays an important role and can have great effect on accuracy and robustness on the VOF 
method (refs. 12 and 13). The main difference between the various surface reconstruction algorithms 
is the assumed interface geometry within cells. Piecewise constant interface calculations, also called 
simple line interface calculation (SLIC), approximate the interface with a straight line aligned with one 
of the mesh coordinates (fig. 4) as opposed to piecewise linear interface calculations (PLIC), which 
allow a nonzero slope for the straight line (fig. 5). 


One can see from figures 4 and 5 that the two VOF methods reconstruct different surface geometries 
from the same volume fraction data. The piecewise linear reconstruction algorithm yields a more accurate 
interface approximation which can lead to a smoother volume fraction distribution in translation and 
rotation problems. This can be important in numerical techniques where surface dynamics (i.e., surface 
tension) is included (like in RIPPLE) and depends on an accurate computation of surface normals, hence 
a “physically smooth” VOF distribution. 


RIPPLE uses a piecewise constant/stair-step VOF advection, called the Hirt and Nichols (HN) 
method. This means the surface is allowed to align with more than one mesh coordinate within a 
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Figure 4. SLIC aligns the surface in a cell with one of the cell coordinates (shaded regions represent 
reconstructed interfaces and numbers shown are volume fractions). 
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A continuous surface line is not obtained 


Figure 5. PLIC allows a nonzero slope for the surface line (shaded regions represent reconstructed 
surfaces and numbers shown are volume fractions). 

cell depending on volume fractions in neighboring cells and the local velocity (fig. 6). For a detailed 
discussion of this advection method see references 5—7. 

For most flows, the HN method gives better results than a SLIC algorithm, but although it is more 
efficient, it is not as accurate as a PLIC VOF advection. In our computations, the surface remains smooth 
and does not undergo severe distortions, where the shortcomings of the HN advection, in comparison to 
the PLIC method, could be crucial. However, it is important to keep the type of surface reconstruction 
in mind when viewing results from RIPPLE computations, especially when surface smoothness, surface 
curvature, and resulting surface dynamics are viewed on a “microscopic” scale of a few cell widths. 

The main advantage of the VOF technique over older surface tracking models like the marker and 
cell (MAC) method or the Lagrangian incompressible (LINC) technique, which use logically connected 
Lagrangian points instead of characteristic marker data, is that the VOF approach can track highly 
distorted and disconnected surfaces. 

Reconstructed surface (stair-step) 



Surface 


Figure 6. The piecewise constant/stair-step algorithm in RIPPLE allows the surface to stair step and 
align with more than one mesh coordinate within a given cell. 
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2.3.2 CSF Model 


The surface-tension-induced pressure jump A p given by equation (8) represents a boundary con- 
dition in the pressure field at the interface. In order to set the correct pressure at the free surface, one 
would have to compute the exact surface location and slope, which would require a time consuming 
reconstruction of a continuous surface. However, a fast and efficient algorithm for reconstructing a 
continuous surface line from VOF data has not been devised yet (ref. 7). An additional problem is that 
the interface line rarely coincides with grid points where boundary conditions such as the pressure jump 
can be easily applied. In addition to a continuous surface reconstruction, an interpolation scheme for 
setting the correct pressure in cells adjacent to surface cells would be required. 

A very different approach to surface forces, in general, is given by the CSF model. It was developed 
by Brackbill, Kothe, and Zemach (ref. 9) and it makes very effective use of VOF data. In the following 
subsection, we will discuss the CSF method briefly to illustrate how it is derived and what the advantages 
are. For a detailed derivation and a proof of the CSF method for surface forces acting normal to the 
interface, see references 9 and 14. 

In the CSF model, a surface tensile force per unit interfacial area (such as surface tension) is 
replaced by a localized volume force acting on fluid elements located in a thin transition region of a 
few cell widths. In the transition, the fluid properties are thought to change continuously from one fluid 
(or void) to another. This means surface forces are no longer applied as discrete boundary conditions 
at a discontinuity (the free surface) but as body forces F v acting everywhere in the transition region. 
Hence the exact surface location is no longer required and F v is applied at grid points located in the 
transition region. 

The basic idea is that one can define a localized body force F v , which volume integral in the limit 
of infinitesimally small transition-region-thickness h is equal to the surface integral of the tensile force 

F s , 

ti>Lv f ' , ^ dv ■ L Ps(x ~* )dS (n) 

where aT s is a position at the free surface. The surface integral is taken over a part of the free surface 
lying within the volume of integration 8V. The orientation of 8V is such that its sides are perpendicular 
to the surface. Figure 7 shows an example using a cylindrical integration volume. In equation (11), 
F s is the classical surface tension force, 

F s (x s ) = crK(x s )n(x s ) (12) 

where h(x s ) is the unit vector normal to the surface. F s (x s ) corresponds to the surface pressure A p 
acting on S: |F 5 (:r s )| = Ap(x s ). 

In order to find the body force F v (x) that satisfies equation (11), first a characteristic function c(x) 
with the following properties is defined (fig. 8) 

ci in fluid 1 

c(x) = < C 2 in fluid 2 (13) 

( c l + C 2 ) /2 at the interface 
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Fluid (F=l) 


Figure 7. Cylindrical volume of integration used for the CSF model. 



Figure 8. Characteristic step function c(x). (Note: x is normal to the surface.) 


From equation (13), one can see that c(£) is a step function and that the interface is specified by 
the discrete value c(x s ) = (ci + c 2 )/ 2 = (c). Hence the interface is an infinitesimally thin line. 

In order to “create” a transition region between the two different fluids (marked by ci and c 2 ),the 
interface is imaginary broadened to a small distance h by convolving c(x) over h with a j B-sphne^ This 
yields a mollified characteristic function c{x) which varies smoothly from ci over (c) to c 2 , as 
in figure 9. The B-spline is defined so that c(x) approaches c(x) as the transiuon region thickn 

goes to zero: lim/ l _o = c (^)- 



Figure 9. Mollified function c(x). 
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With c(x), a bell-shaped distribution function for the volume force F v is derived by taking the 
gradient of c(x) and normalizing with [c] = C2 — c\. Figure 10 shows the distribution along a surface 
normal. The volume force F v is then formulated as 


F v (x ) = aK{x)^Y^- 

l c J 


(14) 


To show that the expression on the right hand side of equation (14) satisfies the identity in the 
integral equation (11), one can apply the following relations. 


Because lim^o c(x) = c(x) (which is a step function) and the normalization of c(x) with (c2— ci), 

lim = H (n(x s ) • (x — x s )) (15) 

h-*0 [cj 

where H is the Heaviside unit function. With equation (15), 

lim Lr aK ( S ) ~Tr’ dV = f° K (&)VH(h{x a ) ■ {x-x s ))dV 
h-*0 J6V [CJ J6V 

Further, 


f aK(x)VH(h(x s ) • (x — x s ))dV = f aK(x)h(x s )6(h(x s ) • (x — x s ))dV 
J6V Jsv 


where S(x) is Dirac’s delta function, which is the derivative of the Heaviside unit function and has the 
property 



6(x)g(x)d(x) = g( 0) 


g(x) is any suitable continuous function and in our case, g(x) = <r/c(x)n(x s ) and <?(0) = cr/c(x s )n(x s ). 
Thus the delta function converts the volume integral into a surface integral. 


/ aK(x)n(x s )6(h(x s ) ■ (x — x s ))dV 
JSV 


/ aK(x s )fi(x s )dS 
J6S 


(16) 


The integrand on the right hand side of equation (16) is the surface force defined in equation (12). 



Figure 10. Vc(x)/[c] along the surface normal; Fj + i, Fj + 2 and Fj + 3 illustrate the discrete values of 
F v applied at grid points i + l,i + 2, i + 3ina numerical calculation. 
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Co mmen ts on F v and k(x)— It is interesting to note that Fy is not a function of the surface 
location x s ; rather, it is a true volume force that depends on x. The discrete boundary value problem 
is replaced with a continuous model. Since c(x) is constant in each Jluid and changes only within 
the transition region, F v becomes zero outside the interface region. F v can be simply added to the 
momentum equation (5) setting F^ = F v . In a numerical calculation, Fy is obtained at grid points that 
are located in the transition region (fig. 10). 


Because nested surfaces of constant c in the transition region represent the interface in the CSF 
model, the normal vector of the interface is obtained by taking the gradient of the mollified characteristic 


function, 

n(x) = Vc(x) 

Substituting the right hand side of equation (17) into equation (7) yields for the curvature 


k ( x ) = — V 


Vc(x) 


(17) 


(IS) 


In RIPPLE’s implementation of the CSF model, the characteristic function c(x) is set equal to the 
VOF function F(x), 


c(x) = F(x) = { 


ci = 0 

> ci,< c 2 

C2 = 1 


in the void, 

in the interface region, 

in the fluid 


since F(x) offers all the properties postulated for c(x). 


233 Wall Adhesion in the CSF Model 

It is obvious from Young’s equation (9) that if the fluid interface forms a contact angle 9 with a 
rigid boundary that is different from the equilibrium contact angle 6 S , a resulting force should be applied 
in order to restore the fluid to the equilibrium contact 6 S . This force can easily be computed within the 
framework of the CSF technique by applying a wall adhesion boundary condition as follows. 

The surface normal vectors at the wall are not computed from equation (17). Instead, the unit 
normal vectors are specified so that they form an angle with the wall normal that is equal to the 
equilibrium contact angle (fig. 11). This is done prior to evaluating the surface curvature (eq. (7)), and 
hence the volume force F v in equation (14) in cells at the wall depends on the difference between 6 
and 9 S . In figure 11, the case for 9 > 9 S is illustrated. One can see that the horizontal component of 
F v is in accordance with Young’s formula, 

O vapor/solid — ^iiquid/solid — ^liquid/vapor COS 9 = <7liquid/vapor( cos 9s ~ COS 6) > 0 


The assumption that the contact angle is constant is a physical approximation in the case of a 
moving contact line (refer to subsection 2.2.3). However, if the difference between the static and 
dynamic contact angle is small (i.e., if the static contact angle and the velocity of the advancing contact 
line is small, like in the current case), this assumption is likely to be a good approximation. 
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Figure 11. Wall adhesion boundary condition in the CSF method; the unit vectors, n bc , at the bottom 
vertices of cells adjacent to the wall are set normal to the equilibrium surface position given by 6 S . 

2.3.4 Application of Tangential Surface Forces 

The driving force of the flow studied in this work (i.e., the flow in the oil patch) is the shear 
force of the air boundary layer acting on the surface of the oil. Since in RIPPLE a solution of the 
flow governing equations is obtained only in the fluid and viscous effects at the surface are neglected, 
a method had to be developed that enabled the application of shear forces at the interface. 

The Viscous Shear Stress Method— It was first suggested to set the shear components of the 
viscous stress tensor in cells containing the surface and empty cells equal to that of the air boundary 
layer, 

{jxy ) oil = (r xy ) air = constant , V cells 0 < VOF <1 (19) 

Figure 12 shows how the shear stress of the air boundary layer is applied in the RIPPLE compu- 
tation. This simple approach is a coarse approximation if the surface deviates from a straight line, has 
nonzero slope, or does not coincide with grid points. It was used to investigate RIPPLE’S applicability to 
our case since the surface was thought to have a slope smaller than 10° at all times. Its main advantage 
is the straightforward implementation of the boundary condition (eq. (19)) into the existing computer 
code. However, it will be shown later that the resulting flow field is nonphysical and the method was 
discarded. 

Extension of the CSF for Tangential Surface Forces- Applying the correct tangential forces at 
the free surface in a numerical scheme is analogous to applying normal forces, which was discussed 
in subsection 2.3.2. Thus a simular approach is used to overcome the restrictions of the viscous shear 
stress (VSS) method. 

Since the boundary layer flow is parallel to the surface it is convenient to define a local coordinate 
system, the x- and y - axes of which lie in the surface tangential plane and the 2 -axis points in the 
outward normal direction. The two-dimensional case is shown in figure 13. 


13 





Figure 12. The application of the shear stress boundary condition. 



Figure 13. The local coordinate system with the principal vectors e x and e 2 is rotated by <p\ n is the 
surface unit normal vector pointing into the fluid. 


The components of the boundary layer shear stress tensor f are given in the local coordinate system 
and can be transformed to the global coordinate system with the tensor transformation, 

Tij = c ik Cji f ki i , j, A;, l = 1 , 2, 3 ( 2 °) 

where c ik and Cji are the transformation coefficients, c ik = cos(Z e t , E k ) and Cji = cos(Z ej, l{). 

The force per unit interfacial area resulting from the shear stress tensor r acting on a surface 
element is then given by 

fV(f s ) = -r(2.) ■ < 21) 


This tensile force can be converted to a localized volume force acting everywhere in the transition 
region utilizing the derivation used in subsection 2.3.2 to obtain the surface tension volume force. 
Setting F(x s ) = F r (x.) in equation (11) yields for the volume shear force 


Frv(x) 


-t{x) ■ 


Vc(x) 

“W" 


( 22 ) 


Figure 14 illustrates how the volume shear forces are distributed in the transition region when the 
surface is oriented horizontally. 

All properties of the volume force F v mentioned in subsection 2.3.2 apply, ofcourse, also to F TV . 
Hence the body force F b in the momentum equation (5) is calculated from F b = F v + F TV . 
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Figure 15. Location of the variables in a grid cell. 
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In the first step of RIPPLE’S projection method, an intermediate velocity field v is obtained from 
jf 1 by dropping the implicit pressure term in equation (23), 


i^ = -V.(W) " + Lv. £ " + r*+L#f (24) 

Using the intermediate velocity field v, the second step becomes an expression for the new velocity 
field u n+1 and the new pressure field p n+1 . 


zl = _ J-v P n+1 

At p n F 


(25) 


Multiplying each side of equation (25) with the Nabla operator and recalling the continuity condition, 
V • v u+1 = 0, one gets a Poisson equation for the pressure. 


V- 



V • v 
At 


(26) 


The solenoidal velocity field v n+l can then be recovered from equation (25). 

It should be noted that the momentum advection terms in equation (24) are estimated using the 
weakly monotonic, second-order upwind van Leer method. The solution of the implicitly discretized 
pressure Poisson equation (26) is obtained through an incomplete Cholesky conjugate gradient technique 
(ref. 16). A detailed presentation of the discretization methods used in RIPPLE is given in references 5 

and 7. 


The explicit treatment of the momentum and VOF advection, the viscous terms, and the body 
forces in equation (23) requires stability time step constraints. Hence for the advection terms a Courant 
condition. 

At < min ^ C 

is applied. The viscous time step constraint is evaluated as 


Ax 


A\u\ 


, C 


Ay 


A|v| 


At < 


(Ax) 2 (Ay)' 


Zv 


(Ax) 2 + (Ay) 2 J 


min 


(27) 


The maximum allowable value for At is also restricted by a requirement which arises from the 
resolution of the propagation of capillary waves on the surface (ref. 7). Thus the time step constraint, 
resulting from the explicit treatment of surface tension forces in the CSF model for a uniform, rectangular 


grid, becomes 


At < 


p(AxAy) 3 / 2 

47TCT 


1/2 


min 


It will be shown in the next section that, for our case, the most limiting time step restriction was 
found to be the viscous time step constraint (eq. (27)). 
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3 COMPUTATIONAL RESULTS 


This section is divided into three subsections. The first presents a brief discussion of initial oil shape 
and grid dimensions. Subsection 3.2 compares the previously discussed techniques for applying shear 
forces at the free surface, namely the VSS and the extended CSF (XCSF) model. In subsection 3.3, 
computational results of the XCSF, the Lax-Wendroff computation, and the linearized solution are 
analyzed and discussed, and the influence of surface tension on the oil flow is illustrated. 


3.1 Geometry and Grid Dimensions 

The initial oil shape is determined by the static contact angle, the fluid volume, and the gravity 
acceleration. Experiments were performed to provide the static contact angle and the dimensions of 
a typical oil drop used in the FISF technique (ref. 4). A contact angle of about 3° was obtained. 
The measured geometric dimensions are shown in figure 16. The complete oil properties are listed in 
Appendix A. 

Since RIPPLE expects an analytical function for the initial fluid distribution (ref. 7), the initial oil 
shape is approximated with a circular arc as shown in figure 16. The resulting drop shape is in close 
accordance with the shape seen in the above mentioned experiments. 

Various grid configurations were tested and the grid sketched in figure 16 was found to give the 
best compromise of resolution and computation time. Hence all computed cases, unless noted otherwise, 
are performed on a rectangular, uniform grid with a grid spacing of 0.025 mm and 0.01 mm in x and 
y directions, respectively. This yields a grid of 322 cells in width and 14 cells in height, including the 
ghost cells surrounding the actual physical domain. 

Due to the very small cell widths (Aar = 0.025 mm, Ay = 0.01 mm) of the grid used, the most 
restricting time step constraint for all computations is the viscous time step constraint (eq. (27)). With 
the given grid spacing, it becomes At = 5.747 • 10 -7 s. 



Figure 16. Initial oil distribution and grid dimensions. 
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3.2 Comparison of the VSS Method and the XCSF Model 

To compare the two methods for applying tangential surface forces, the flow discussed in subsec- 
tion 2.1 is computed. Two different cases are presented in which the surface tension coefficient a is 
varied. The first case is the “standard” case where a = 0.0208 N/m was set (the a measured by e 
oil manufacturer). In the second case, a was set to a value three times higher. The value for the shear 
stress t = 2.394 N/m 2 , as measured using the FISF in a channel flow facility (ref. 4), was used. 

When viewing the following figures, it is important to keep in mind that the y-axis is scaled by a 
factor of roughly 40 compared to the x-axis. Hence the oil layer is much shallower than suggeste m 
figures 17—19 (and in figures 21—23 in subsection 3.2.2). 


3.2.1 The “Standard” Case 


In figures 17-19, the evolution of the drop shape computed by the XCSF and the VSS^is shown 
at three different times. The initial shape as well as the applied shear stress, r = 2.394 N/m , and t e 
surface tension coefficient a are the same for both computation methods. From figures 17-19, one can 
see that the drop shape can be divided into two main parts: the front part with the leading edge an t e 
downstream part containing the trailing edge. 


The Leading Edge- At early times, the leading edge of the drop is deformed to a wedge, with the 
leading edge angle that decreases and a wedge length that increases with time. In the wedge area, the 
XCSF results show a slightly smoother surface than the VSS computation, where the surface line tends 

to “stair step” (figs. 18 and 19). 


It is interesting to see that the “step height” is of the same order of magnitude as Ay. The “stair 
stepping” of the surface line is amplified when the interface approaches a horizontal orientation. 

This difference in surface quality between the VSS and the XCSF method results from the refined 
application of surface shear stress forces in the XSCF. The VSS method does not take surface slope and 
surface sub-cell location into account as opposed to the XCSF method. In the VSS technique, the shear 
stress is set to r xy = 2.394 N/m 2 in cells containing the surface and empty cells (i.e., in cells where 
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Figure 19. Drop shape evolution at 2.4 s, a = 0.0208 N/m. 


0 < VOF < 1). Since in RIPPLE the off-diagonal components of the viscous stress tensor are computed 
at cell vertices, the expression for the distance d between the location of the surface boundary condition 
T xy = 2.394 N/m 2 on the grid and the surface line becomes 0 < d < Ay. This means that although the 
distance of the application point of the surface boundary condition to the actual surface within a cell is 
a transient function, the value of the VSS stress boundary condition is held constant. Furthermore, the 
application point is a noncontinuous function; it “jumps” by Ay as soon as the fluid surface “swaps” 
into the next cell row. 

Hence the stair-stepped character of the wedge computed by the VSS technique is thought not 
to be of physical origin and is therefore, obviously, an undesired behavior (with respect to the FISF 
technique). In the FISF, the thinning rate of the oil wedge is measured at a fixed z-location and it 
is assumed to be inversely related to time (see equation (3) and reference 4). However, the wedge 
thinning rate obtained from the stair-stepped surface line in the VSS computation does not show a linear 
behavior: each time a surface line step passes through the observation point, the height decreases more 
rapidly than in regions away from the step area. 

The Trailing Edge- A significant difference in the two methods can be seen near the trailing edge 
of the drop, where the curvature (and hence surface tension) is higher than at the leading edge. At 
0.4 s, the difference is still hardly visible, but figures 18 and 19 show that it becomes more pronounced 
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at later times: the downstream contact line of the oil drop in the XCSF solution moves faster than the 
one computed using the VSS method. 

Since the actual contact angle at the downstream edge of the drop is greater than the specified static 
contact angle. Young’s formula (eq. (9)) yields a force component that tries to advance the contact line 
downstream and therefore works in addition to the applied surface shear stress. It is surprising that the 
difference in the velocity of the advancing contact line is so high, since m both methods the contact 
angle is almost of the same value, and hence nearly the same accelerating force component acts on the 

trailing edge contact line. 

The lower velocity of the downstream contact line in the VSS computation in comparison with 
the XCSF, in combination with similar thinning rates of the wedge in the front part of the drop, results 
in a greater maximum oil height and a more highly curved surface in the VSS result. It seems that 
the energy, which is introduced to the system through the external air flow, is to a great extent used to 
deform the drop rather than accelerate it. In the VSS method, more energy is converted into surface 
tension energy than into kinetic energy in comparison with the XCSF results. 

Considering the units of F v (x) (surface tension volume force) as energy per unit volume [J/m 3 ] 
(eq. (14)), a comparison of the total amount of surface tension volume forces (the sum of surface tension 
volume forces in all cells) acting on the oil illustrates the difference in surface tension energy. One 
important point to recall is that the total amount of fluid specified at the beginning of the computation 
is the same in both techniques. 

Figure 20 shows the comparison of F v (x) as a function of time. It can be seen that both functions 
start out at the same point, since the initial shape is the same for both computation methods. However, 
the total amount of surface tension volume forces, and hence surface tension energy, is c ^ ea ^ y § rea ^ er 
in the VSS case at all computed times. This leads to the conjecture that a greater retarding force must 
act on the fluid in the VSS method, which prevents the trailing edge of the oil from moving downstream 
at as high a rate as the XCSF drop and causes the high surface deformation. 
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In order to verify that the amount of surface energy (i.e., surface tension volume force) due to 
the high surface deformation in the VSS technique leads to a lack of kinetic energy in the fluid, a 
kinetic energy balance was performed: the square of the magnitude of the velocity vector in each cell 
was multiplied with the corresponding value of the VOF function F. (Note: F reflects the fractional 
volume, and hence fractional fluid mass, in a single cell.) These “scaled” velocity magnitudes were 
then summed over all N grid cells, 

N 

e kin = ' M i (28) 

i=l 

which yields the specific kinetic energy in the fluid. In table 1 the numerical values for e^ n are listed 
for t = 2.4 s; for other times, the difference in the value for e t t-i n is similar and will be not listed. 


Table 1. Values for e kin at t = 2.4 s. 


Time 

^ kin VSS in 

mm 2 /s 2 


e kin XCSF in 

mm 2 /s 2 


2.4 s 

1034.65 

2696.53 


3.2.2 The “High Sigma” Case 

A lower contact line velocity yields a lower capillary number ( Ca = uji/ cr) which in turn means 
that interfacial forces play a more important role in the VSS technique. Hence the second case with 
cr = 0.0624 N/m (figs. 21—23) was computed to investigate the impact of surface tension on the velocity 
difference at the contact line and finally to determine which shear force method yields the more physical 
flow. 


The results are basically the same as in the “standard” case described in subsection 3.2.1 with the 
difference being that the trailing edge contact line in the VSS computation does not move at all. This 
seems to be a paradox, since the contact angle is again the same in both computations, and hence the 
wall adhesion force should move the contact line of the oil drop downstream in the VSS computation 
as well. This is an illustration of the physical deficiencies of the VSS implementation. 



Figure 21. Drop shape evolution at 0.4 s; cr = 0.0624 N/m. 
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In figures 24 and 25 the velocity profiles of the two techniques, in the area of greatest oil height, 
are shown. One can see that the velocity field of the oil flow in the XCSF computation loo s very 
simular to a Couette flow, which was expected. 

An unexpected result of the VSS method is that the oil in the vicinity of the surface flows in the 
same direction as the external air flow, but in the oil the fluid moves upstream against the air flow 



red: Surface Line 
black: Velocity Vector* 


Figure 24. Velocity vectors XCSF computation, o = 0.0624 N/m, t - 3.4 s. 
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red: Surface Line 
black: Velocity Vectors 


Figure 25. Velocity vectors VSS computation, a = 0.0624 N/m, t = 3.4 s. 


direction. This means that the oil flow circulates inside the drop instead of moving downstream with 
the air flow. The explanation for this unexpected flow behavior can be found through analyzing the 
VSS method by considering the action of the surface tension forces. 


Setting the shear stress at the surface to a constant value (like in the VSS method) forces the fluid 
into a motion with velocity gradients. Or one can say that the boundary condition is satisfied when the 
flow field possesses velocity gradients at the surface that fulfill the expression 


N 


du dv 


TXy ~ 2 394 77 ? ~ + & 



Since the initial drop shape is the state of minimum surface tension energy (fig. 20), any deviation 
from this state induces surface forces that try to restore the initial surface position and hinder the 
acceleration of surface fluid elements. Hence the velocity gradient resulting from the shear stress 
boundary condition is established not only by accelerating fluid elements at the surface in the downstream 
direction, which was the desired effect, but also by moving fluid particles inside the flow in the upstream 
direction. This flow field satisfies all given boundary conditions but certainly does not reflect the physical 
facts as observed in the Fluid Mechanics Laboratory at Ames. 


3.23 XCSF Vectors 


Figure 26 shows the oil drop surface and the resulting volume shear force vectors in the region 
of highest curvature in the “standard” case. (It should be mentioned that the x- and y-axes are scaled 
equally in this plot.) The width of the transition region and the alignment of the XCSF vectors tangential 
to the surface can be seen very well and seem to be physically correct. 

An additional comment should be made at this point about the computation of F TV (x). In RIPPLE, 
equation (22) is implemented with an additional term, 

F tv (x) = -t(x)VF(x) ■ 2F(x) (29) 


23 



Red: Surface Line (at VOF-O.5) 
Black: Volume Shear Force Vectors 


Figure 26. Volume shear forces. 


The multiplication by 2 F{x) is permissible, since it does not affect the validity of 

lim [ F TV (x)dV = [ F T (x- s )dS 
h^OJSV J6S 

which is the equivalent expression to equation (11) for the volume shear force. This is obvious bec ^*se 
through the definition of the mollified function (i.e., F(x) in RIPPLE), the limit expression (eq. (15)) 

becomes p , p 

lim [F(x) ■ 2F(x)] = 2 1 - lim F(x) = H{n{x s ) ■ (x - x s )) 


The additional factor in equation (29) biases the volume force more into the fluid than suggested in 
figure 14. Figure 26 shows how the shear force vectors are more pronounced in the part of the transition 

region that lies inside the fluid. 

The biasing yields the advantage that accelerations due to the shear forces are proportional to VF 
rather than F itself. F VT becomes then a true body force. (A similar procedure for interfacial volume 
forces due to surface tension is discussed in references 7 and 9.) 


3.2.4 Computational Conclusions 

In an effort to extend RIPPLE’S capabilities to model surface shear forces, two different approaches 
were implemented and evaluated. 

One is the VSS method, which simply sets the viscous shear stress in surface cells to the desired 
value. The second approach, which is a more refined stress application technique, is an extension ot 
the CSF model that was derived and implemented. 

Both methods were tested with the “standard case,” which reflects the physical situation of the flow 
channel facility where the 3D-FISF technique was developed (ref. 4). 


It was found that the XCSF model gives more physical results. In particular, the circulating veloci y 
field inside the drop computed by the VSS method is thought to be unphysical. The velocity of the 
contact line at the trailing edge of the oil, resulting from the VSS computations, does not agree wi 


the theory of advancing contact lines. 
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Thus the investigation of surface tension and wall adhesion effects on the FISF technique in the 
following subsection is carried out using the XCSF model. 

33 Comparison with Previous Studies; Surface Tension and Wall Adhesion 

Effects 

In this subsection, we will first present h(t) graphs obtained from the Lax-Wendroff algorithm and 
the linearized equation (3), which is used in the FISF technique. 

Subsequently, the impact of wall adhesion is investigated through the variation of the static contact 
angle. The effects of surface tension are illustrated through computational results obtained using different 
values of the surface tension coefficient. 

3.3.1 Lax-Wendroff Results and the Linearized Equation 

Figure 27 shows the h(t ) plots of the linear solution and Lax-Wendroff computation at a distance 
x = 1-5 mm from the leading edge. One can see that the plots lie on top of each other after about 
10 s. At this time, the fluid height is approximately 3 • 10~ 3 mm. In order to resolve the flow in an oil 
drop of comparable height using RIPPLE, a much finer grid than described in subsection 3.1 would be 
necessary. Since the limiting time step in all RIPPLE computations was determined by the viscous time 
step constraint, the maximum allowable time step decreases by the power of 2 with a grid refinement 
(eq. (27)). In addition, the number of grid points would increase quadratically. 

Thus a grid that resolves the oil flow in a 3 ■ 10 -3 mm high drop would yield a tremendous increase 
of computation time. Currently, the CPU time for a typical computation of the drop evolution up to 
2.4 s amounts to 108 CPU hours on an SGI R10000 machine. 



Figure 27. Lax-Wendroff and linearized equation, x = 1.5 mm. 
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Since surface tension forces are localized body forces, it is assumed that surface tension effects are 
limited to regions of high curvature. From the drop evolution figures in the preceding section, one can 
see that the area of high curvature is located near the trailing edge, which moves downstream y 
from the region of interest as far as the FISF measurement technique is concerned. Hence the effects ot 
surface tension on h(t) at x = 1.5 mm are expected to decrease with time and the time frame of 2.4 
is thought to be sufficient for the investigations. 

An additional note must be made about the “high sigma” computation (<r = °-062 4 N/m). Although 
the resulting oil height at x = 1.5 mm from the leading edge is lower than one cell width the oil flow 
was computed up to 3.4 s. This was done in order to compare the global behavior of the oil drop flows 
computed with the XCSF and VSS methods only (refer to subsection 3.2). The behavior of die flow m 
the wedge area, and in particular h(t) for x = 1.5 mm, is not expected to be resolved sufficiently at 
times greater than approximately 2.4 s. Hence it is not referenced m the following discussion. 

33.2 The Influence of Wall Adhesion 

In order to determine the effect of wall adhesion on the oil flow, two cases were computed with 
RIPPLE, where the contact angle was set to 9$ = 0.1° and 9 S — 6°, respectively. (<r was speci e o 
0.0208 N/m in both cases.) 

The relatively low velocity of the contact line suggests that the dynamic contact angle does not 
differ significantly from the static contact angle. Thus the given values for 6 S are thought to be extrema 
which are well suited to answer the question of the influence of wall adhesion on the drop motion, and 

hence on the FISF technique. 

Figure 28 shows the results of the two cases. The computed drop shapes are almost identical at all 
plotted times. In particular, the surface lines coincide in the wedge area, which is the mterestmg jT“ 
of the oil drop as far as the FISF technique. Even in the vicinity of the two contact lines, no sigmfican 

differences are visible. 

This result is not surprising, since it follows from Young’s formula (eq. (3)) that the resulting 
horizontal force acting on the contact line in the case of nonequilibrium is proportional to cos0„ which 
varies between 0.999 and 0.995 for the specified angles. This means that the difference between 
contact line forces in the two cases is marginal and has no significant effect on the oil ow. 
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Two conclusions can be drawn from the results shown in figure 28: (1) wall adhesion forces 
do not affect the oil drop motion for a wide range of 6, and (2) the simplification made in RIPPLE, 
that the contact angle at the moving contact line is constant and equal to the dynamic one, is a valid 
approximation for the flow problem studied in this work. 

3.33 Surface Tension Effects 

RIPPLE Compared with Lax-Wendroff- Figure 29 displays the drop evolution in the wedge 
area obtained from the Lax-Wendroff computation and RIPPLE. The spatial restriction to the front part 
of the drop arises from assumptions implicit in the implementation of the Lax-Wendroff algorithm. It 
should be noted that the ^-direction grid spacing used for these two computations is the same. 

One can see that the oil quickly forms a wedge in both computations. The slope of the surface line 
computed with the Lax-Wendroff technique decreases faster at early times than the one obtained from 
RIPPLE. This behavior is not surprising, since surface tension was identified to be a retarding force that 
tries to establish the initial shape and is omitted entirely in the Lax-Wendroff calculation. 

Variation of cr in RIPPLE- The retarding effect of surface tension on the oil motion is also 
visible in figures 30-32, where the drop evolutions from RIPPLE computations for a = 0.0208 N/m 
and a = 0.0624 N/m are shown. 

The trailing edge moves slower in the downstream direction and the surface line is less curved 
for the case with a = 0.0624 N/m. The point where the drop shape deviates from the wedge form is 
located more upstream than in the “standard” case. 

In figure 33, surface tension volume forces are plotted as functions of time for the two cases. It 
is interesting to see that although the surface line is less curved in the “high sigma” case, more energy 
per unit volume is used for the surface deformation. This explains the slower downstream motion of 
the oil drop in the high sigma case: since the flow energy introduced to the oil through the external air 
flow is the same at t = 0 in both computations, the higher surface tension energy causes a reduction of 
kinetic energy. 



Figure 29. Drop evolution (wedge area), cr = 0.0208 N/m, 6 S = 3°. 
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It looks like surface tension not only affects the surface line curvature but also the overall fluid 
flow. This seems to be in contradiction with the statement made previously that surface tension forces 
are localized body forces with local effect on the fluid flow. 

Looking at the global drop shape, all computations show a similar flow behavior: the back part of 
the drop surface is more or less curved while the front part becomes a wedge with a straight surface 
line After a certain period of time, all computed drops show quite a uniform downstream motion 
which means that the thinning rate of the wedge, h(t) (which is measured in the FISF technique) does 
not differ fundamentally when surface tension is included in the calculation. It is bought that once Ae 
uniform state of motion is reached, the change of h(t) by including surface tension canbe extrapolated 
in time in order to estimate the total effect of surface tension on the FISF technique. This is presented 

in the following subsection. 


28 





The Influence on h(t ) and the FISF Technique— The influence of surface tension on the wedge 
thinning rate is illustrated in figure 34, where the RIPPLE results for a = 0.020 8N/m and o = 
0.0624 N/m as well as the Lax-Wendroff results are included. 

The h(t) graphs lie on top of each other until about 0.2 s. After this time, the surface deformation 
is high enough so that the resulting surface tension force restrains the drop motion, and hence retards 
the evolution of the wedge. One can see that increasing a in the RIPPLE computations yields a greater 
wedge height at the same location and time. 

At about the 0.6 s point, the various solutions for h(t) approach straight lines in the logarithmic 
scaled plot of figure 34. After this time, RIPPLE’s h(t) graphs remain parallel to the Lax-Wendroff 



Figure 34. h(t) at x = 1.5 mm. 
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graph, which approaches the linear equation, 

.... Moil* _ constant 
T wa iit t 


as shown in figure 27. 

Functions of this form (having different y-values but the same slope in a logarithmic scale) differ 
only by a constant factor. Thus an analytic expression for h(t) obtained from RIPPLE can be derived 
when “scaling” the Lax-Wendroff solution with an appropriate factor /, 


h(t)R!PPLE — / ’ h(f)Lax-WendrofT / 


constant 


Vi > 1.5s 


With this equation, the relative height difference between RIPPLE’S results (including surface 
tension) and the Lax-Wendroff computation, 


drel 


^(t) RIPPLE — h(t)Lax-Wendroff 
fi(t)RlPPLE 


/-I 

/ 


becomes a constant expression rather than a function of time. 

To determine /, the heights at the latest possible time are used, in order to be as close as possible 
to the point where the Lax-Wendroff computation coincides with the linear expression (fig. 27). Since 
the computation with cr = 0.0624 N/m is a numerical experiment not reflecting the physical situation, 
/ was calculated only for <7 = 0.0208 N/m. The results are summarized in table 2: 


Table 2. Values for / and d re i at t = 2.4 s. 


Time 

^Lax-Wendroff 

Cripple 

f 

d re i 

2.4 s 

1.07511 • 10~ 2 mm 

1.15685 • 10 -2 mm 

1.0707 

7.07% 


Comments on /- RIPPLE’s results show that the specified shear stress value yields a 7% smaller 
oil thinning rate when surface tension is included than predicted by * e Lax-Wen^off computation and 
the linearized solution. In order to obtain the same oil thinning rate with RIPPLE, one would have to 
specify a 7% higher shear stress value. This means that the shear stress level measured with te 
technique (which uses the linear expression to relate shear stress to oil height) is underpredicted by 7%. 
However, the value for / listed in table 2 is based on data that are taken at a very early time ( . s) 
compared with the typical time in the FISF technique (e.g., 100 s). In addition, it was assumed tha .the 
linear oil thinning behavior in RIPPLE after 1.5 s persists, which is not a certainty since the oil flow 
is governed by nonlinear equations. There is limited evidence measured in the laboratory that indicates 
that the linear oil thinning behavior does not persist at later times. If this proves to be the case then t e 
actual underprediction of r wall may be much less. Thus 7% underprediction for r wall by the FISF is a 
very rough estimate, which requires more investigation. An additional uncertainty arises from 
that in RIPPLE only the oil flow is modeled and the effect of the surrounding air flow is reduced to the 
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assumption that shear forces at the interface are the only way the air flow acts on the oil. Any further 
interaction between the oil flow and the air boundary layer (e.g., pressure gradients) is omitted. Thus 
the value for / should be interpreted as a more qualitative statement rather than an absolute “correction 
factor” for the FTSF technique. 

It is also expected that /, and hence the influence of surface tension, is less for larger volumes of 
oil. By enlarging the oil drop, surface line curvature decreases while the area where shear forces are 
applied increases. This means that surface tension forces are weakened and the viscous forces at the 
interface are amplified. 

3.3.4 Influence of Grid Resolution 

The purpose of the last subsection of this section is to verify the adequacy of the chosen grid 
resolution. For this purpose, the previously described “standard” case was computed on a grid, where 

Ax = 0.05 mm and Ay = 0.02 mm. This corresponds to a doubling of the grid spacing in both 
directions. 

In figure 35, h(t) is plotted for the two different grids. The two graphs lie on top of each other 
until the oil height in the coarse grid case is significantly less than one cell height. At that time (i.e., 
approximately 2 s), the oil height begins to decrease more rapidly in the computation with the coarse 
grid. The fact that the oil thinning rate in the case with the doubled grid spacing is almost the same as 
on the finer grid (until about 2 s) shows that the grid in the previously discussed computations is fine 
enough for our purposes. 



Figure 35. The influence of grid resolution on h(t), x = 1.5 mm. 
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4 CONCLUSIONS AND RECOMMENDATIONS 

(1) A computational study was performed to investigate the effects of surface tension and wall 
adhesion on the accuracy of the FISF technique for the measurement of skin friction. 

The computer code RIPPLE, which solves the two-dimensional, incompressible, transient Navier-Stokes 
equations for free surface flows, was used to compute the flow of an oil line placed on a flat pla e in 
an air boundary layer. RIPPLE was modified in order to model shear forces acting on the free surface^ 
In this context, two different approaches for applying shear stresses at the interface were deve ope 
and discussed. One is the viscous shear stress model, where the value for the skm friction measured 
with the FISF technique is set explicitly in grid cells containing the interface. The second approac 1 
an extension of the continuum surface force model (CSF) for tangential forces, where skin faction is 
converted into a localized body force and simply added to the momentum equations. 

The two models were compared and it was shown that the extended continuum surface force model is 
ideally suited for the application of tangential surface forces (i.e., shear forces). It was used for the 
subsequent study of surface tension and wall adhesion effects on the oil drop motion an en 
accuracy of the FISF technique. 

(2) The surface tension coefficient a was varied to illustrate the effect on the oil drop ^ 

deformation. Different static contact angles were specified to investigate the influence of wall adhesion. 
In all case studies, special focus was made on the thinning rate of the leading edge o e 01 p ’ 
since this is the drop area where the FISF technique is applied. In addition, RIPPLE s result for the oi 
thinning rate was compared with previous studies discussed in reference 4. 

(3) It was demonstrated that the variation of the contact angle does not affect the flow of the oil 
drop, and hence wall adhesion has no impact on the accuracy of the FISF technique. 

The variation of the surface tension coefficient and the comparison with Zilliac s work (ref. 4) showed 
that surface tension has a retarding effect on the oil thinning rate, which may lead to an underprediction 
of the skin friction measured with the FISF technique. This effect can be reduced y en arging P 

volume, which yields a lower surface curvature and a higher sum of driving tangential shear forces 

the interface. 

(4) Finally, some recommendations for improving the results presented in this paper and reducing 
the uncertainty of the FISF accuracy estimation are now made in terms of future work: 

• Establish an implicit treatment of the viscous term in the momentum equation in order to 
overcome the restricting viscous time step constraint and to be able to compute further in 

time on a even finer grid. 

• Extend RIPPLE to two fluids in order to model the entire physics involved in the flow. 

• Include a turbulence model for modeling the turbulent air boundary layer. 

• Extend RIPPLE to three dimensions to validate the 3D-FISF technique. 


Some of these items, such as the capability of modeling two fluids, are currently under development. 
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APPENDIX A 
OIL PROPERTIES 


A.l Dow Corning Silicone Fluid 

Name 200 Fluid 

Manufacturer Dow Coming Corporation, USA 

Density at 25°C 9600 kg/m 3 

Surface tension at 25°C 0.0208 N/m 

Kinematic viscosity at 25°C 50 Cs 


A.2 Description 

The 50 Cs 200 Fluid from Dow Coming is a medium viscosity poly-dimethyl-sil-oxane polymer. 
It has the following typical chemical composition: 

(CH 3 )3SiO[SiO(CH3) 2 ] n Si(CH 3 )3 
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APPENDIX B 
RIPPLE 

The original version of RIPPLE is available through the National Energy Software Center, 
9700 South Cass Avenue, Argonne, IL 60439. 

B.l A Typical RIPPLE Input Deck 


Oil drop 
$numparam 

conserve= . false. , alpha =2 . 0 , 
twfin=2.6, prtdt=500.2, pltdt=0.2, 
delt=0 . 0000001 , dtmax=0 . 1, 
nrestart=-l . 0, 
dmpdt=0 . 2 , 

kb=2 , kt=2 , kl=2 , kr=2 , 
autot=l . 0D0 , 
sym= . true . , 
gfnctn=. true. , 

Send 

$f Idparam 
icyl=0, isurfl0=l, 
gy=-9799 . 3 , 

cangler=3 . 0 , canglel=3 . 0 , 
canglet=3 . 0 , cangleb=3 . 0 , 
xnu=50 . 0 , rhof=0. 00096, sigma=20.8, 
psat=101300 . 0 , 

Send 

$mesh 

nkx=l , xl=0. 0,8.0, xc=4.0, 
nxl=160, nxr=160 , dxmn=0.025, 
nky=l , yl=0. 0,0.12, yc=0.01 f 
nyl=li nyr=ll, dymn=0.01, 

Send 

Sshear 

shear t = 2.394, 

Send 

Shovertlst 
xloc=l . 7d0 , 

Send 

Stimestepadj 
tmf ac=l . 5d0 , 

Send 

Sobstcl 

nobs=0. 

Send 

$f reesurf 
nfrsrf=l, 

fal ( 1) =5 . 4 , fa2(l)=-l, 

fbl ( 1 ) =-95 . 4056834388, fb2 ( 1) =-l , 

fcl(l)=-l-04, 

ifh(l)=l, 


34 



$end 

$graphics 

plots= . true . , dump= . false . , 
Send 


This input deck was used to compute the previously discussed “standard” case. For a detailed 
discussion of the input data, see reference 7. 

B.2 Source Code of Main Subroutines for the Implementation of the XCSF Model 

Two subroutines are listed. The first one serves the purpose of computing the fluid accelerations 
due to surface shear forces. The second one computes the surface shear forces using the XCSF method. 
Many other subroutines were modified or added to the original version of the code but they will be 
omitted in this paper for reasons of space. 


subroutine shrfrce 
c 

c 

c Purpose - 

c compute fluid acceleration due to shear forces 


name 


c SRFFRCE is called by - 
c 

c name 

c RIPPLE 

c 


c SRFFRCE calls the following subroutines and functions - 
c 


BC ripple 


c 

C############### ############################################### 

implicit real *8 (a-h,o-z) 

# include "32bit.h'* 

c################################# ############################# 

c 

c############ 

# inc lude *' comdk2 . h H 
#include "iccgdkl.h" 

c############ 

c 

data tiny /1.0d-25/, zero /O.OdO/ 
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c 


do 20 j =2 , jml 
do 20 i=2, iml 

ij= ( j-1) *imax+i 
ipjsij+1 
ijp=i j+imax 

rhoi j = f ( i j ) *rhof 
rhoipj=f (ipj ) *rhof 
rhoijp=f (ijp) *rhof 

rhorc=(delx(i+l) *rhoij+delx(i) *rhoipj ) / (delx(i) +delx(i+l) ) 
rhotc= (dely ( j + 1) *rhoi j+dely ( j ) *rhoijp) / (dely{ j ) +dely( j+1) ) 
rhox = r ho r c + 1 i ny 
rhoy=rhotc+tiny 

f tauxi j=f taux l i j ) 
f tauyi j=f tauy ( i j ) 
ftauxipj=ftaux(ipj ) 
f tauyi jp=f tauy (ijp) 
c 

c compute the volume shear forces x-component at the right cell face 

f xrc= ( delx ( i ) * f tauxip j *delx ( i+1 ) * f tauxi j ) / ( delx ( i ) -delx ( i+ 1 ) ) 
if (gfnctn) fxrc=2 . *rhox* fxrc/ rhof 
c 

C.... compute the volume shear forces y-component at the top cell face 
c 

fytc= (dely ( j ) * f tauyi jp+dely ( j+1 ) ‘ftauyij ) / (dely (j ) +dely{ j+1) ) 
if (gfnctn) fytc=2 . *rhoy* fytc/rhof 


if (ar(ij) .lt.l.0d-06) go to 10 
u { i j ) =u ( i j ) +delt*fxrc/rhox 

10 u(i j } =cvrogt (zero ,u(ij ) , (rhorc .It. frctn*rhof ) . or. 

& (ar (i j ) . It . 1 . 0d-06) ) 

if (at(ij) .lt.l.0d-06) go to 15 
v(ij )=v(ij ) + de 1 1 * f y tc / rhoy 

15 v ( i j ) =cvmgt ( zero , v (i j ) , (rhotc. It . frctn*rhof ) .or. 

& (at (i j ) .It . 1 . 0 d- 06 ) ) 

20 continue 


c 

c 

c 

c 


call be 


return 

end 


36 



subroutine tension 


c 

c 

c Purpose - 

c Compute the volume forces due to surface tension 
c 

c TENSION is called by - 
c 

c name name name name name name 

C RIPPLE SRFEMIN 

c 

c 

c TENSION calls the following subroutines and functions - 
c 

c name source name source 

c 

c BDYCELL ripple CRVATURE ripple 

c ISMAX cftmath MOLLIFY ripple 

c 

c 

c############################################ ################## 
implicit real*8 (a-h,o-z) 

# include "32bit.h" 

c 

C############ 

# include " comdk2 . h“ 

#include " iccgdkl . h" 

c############ 


name source 

SETARRY ripple 


c 

dimens ion r o ( nxy ) 

dimension vnormx ( nxy ) , vnormy (nxy) 
c 


data tiny /l.Od-25/, zero /O.OdO/ 



call setarry 
call setarry 
call setarry 
call setarry 
call setarry 
call setarry 
call setarry 
call setarry 


( ro, 0 . OdO , nxy) 
(gradrox, 0 . OdO , nxy) 
(gradroy , 0 . OdO , nxy) 
{ tensx, 0 . OdO, nxy) 

{ tensy, 0 . OdO, nxy) 

{ fsv, 0 . OdO , nxy) 

( f tilde, 0 . OdO , nxy) 
{kappa, 0 . OdO, nxy) 


c 



, . . Compute cell densities and copy the 
raw color function, f(ij), into the 
smoothed color, ftilde(ij) 

do 50 j=2,jml 
do 50 i=2 , iml 
ij= ( j-1) *imax+i 
ro(ij)=f (ij)*rhof 
ftilde (i j ) =ro (ij ) 

50 continue 

Reflect ro(ij) and ftilde(ij) in ghost cells 


call bdycellUml. jml, 1,1, 1,1,0. OdO, 0. OdO, 0. OdO, O.ODO, pbc.ro) 

call bdycelltiml, jml, 1,1, 1,1,0. OdO, 0 . OdO, 0 . OdO , O.ODO, pbc, ftilde) 


Optionally smooth the color function 


if (smooth) 

St call mollify ( 2 , iml , 2 , jml , 1 , imax, nsmooth, ro, ftilde, af ) 


Compute vertex- centered normal vectors 

do 35 j=2 , jml+1 
do 35 i=2 , iml+1 
ij= ( j-D *imax+i 
imjm=i j -1-imax 
ijm=i j -imax 
imj = i j -1 
im=i-l 
jm= j -1 

deltay=0 . 5* (dely l j ) +dely ( jm) } 
deltax=0 . 5* (delx (i) +delx ( im) ) 

rhoc= (delx ( i ) *ro ( imj ) +delx ( im) *ro ( i j ) ) / ( delx ( i ) +delx ( im) ) 
trhot= (delx (i) ’ftilde (imj ) +delx(im) *f tilde (ij ) ) / 
k ( delx ( i ) +delx ( im) ) 

rhob= (delx(i) *ro(imjm) +delx(im) *ro(ijm) ) / (delx(i) +delx< im) ) 
trhob= (delx ( i) * ftilde (imjm) +delx(im) *f tilde (ijm) ) / 

& (delx (i) +delx(im) ) 

rhol= (dely ( j ) *ro (imjm) +dely ( jm) *ro (imj ) ) / (dely { j ) +dely ( jm) ) 
trhol= (dely ( j ) *f tilde {imjm} +dely ( jm) * ftilde (imj ) ) / 

St ( dely { j ) +dely ( jm) ) 

rhor= (dely ( j ) *ro (ijm) +dely ( jm) *ro(ij) ) / (dely ( j ) +dely ( jm) ) 
trhor= (dely ( j ) * ftilde (ijm) -t-dely (jm) * ftilde (ij ) ) / 

& (dely ( j ) +dely ( jm) ) 

gradrox { i j ) = ( rhor-rhol ) / del tax 
gradroy (ij ) = (rhot-rhob) /deltay 


c.... store the non-unit surface normal vectors 
c 

vnormx (i j ) =gradrox (i j ) 
vnormy ( i j ) =gradroy ( i j ) 
c 

gradf tx ( i j ) = ( trhor-trhol ) /deltax 
gradf ty ( i j ) = ( trhot-trhob) /deltay 
35 continue 
c 

c.... Impose wall adhesion via an equilibrium 
c contact angle in all obstacle cells 
c 

do 2350 j=2,jml 
do 2350 i=2,iml 
ij= ( j -1) *imax+i 

if (ijobs (i j) .eq.O) go to 2350 
if (ijobs (ij) .eq.l) thetaeq=cangleb 
if (i jobs ( i j ) . eq. 2 ) thetaeq=canglet 
if (ijobs (ij ) .eq.3 ) thetaeq=canglel 
if (ijobs (i j ) .eq. 4) thetaeq=cangler 
c 

grdmgro=sqrt (gradrox (i j ) * *2 + gradroy ( i j ) **2 ) + tiny 
grdmgft=sqrt (gradf tx ( ij ) **2 + gradf ty ( ij ) **2 ) + tiny 
tx= -gradnwy { i j ) 
ty=gradnwx ( i j ) 

csthtar= ( tx*gradrox ( i j ) +ty*gradroy { i j > } /grdmgro 
csthtaf= (tx*gradftx (ij } +ty*gradfty {ij > ) /grdmgft 
txr=cvmgt ( -tx, tx, csthtar . It . zero) 
tyr=cvmgt ( -ty, ty, csthtar . It . zero) 
txf =cvmgt ( -tx, tx, csthtaf . It . zero) 
tyf =cvmgt ( -ty , ty , csthtaf . It . zero) 
c 

grdtwrx= grdmgro * txr 
grdtwry= grdmgro* tyr 
grdnwrx=grdmgro*gradnwx (i j ) 
gr dnwry=grdmgro* gradnwy (ij ) 
grdtwfx=grdmgf t* txf 
grdtwf y=grdmgf t * tyf 
grdnwf x=grdmg ft* gradnwx (ij ) 
grdnwf y=grdmg ft* gradnwy { i j ) 
c 

gradrox ( i j ) =grdnwrx*cos ( thetaeq) +grdtwrx*sin ( thetaeq) 
gradroy ( i j ) =grdnwry*cos { thetaeq) +grdtwry*sin ( thetaeq) 
gradf tx ( i j ) =grdnwfx*cos { thetaeq) +grdtwfx*sin { thetaeq) 
gradf ty ( i j ) =grdnwfy*cos ( thetaeq) +grdtwfy*sin ( thetaeq) 
c 

2350 continue 
c 

c.... Compute the mean curvature, kappa, using 
c the smoothed color function, ftilde(ij) 

c 

call crvature(2, iml, 2, jml, 1, imax, f tilde, gradf tx, gradf ty, 

& delx, dely, r, ri, kappa, tiny) 

c 

c 
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do 110 j = 2 , jml 
do 100 i=2,iml 
ij = ( *imax+i 

ipj=ij + l 
imj =i j -1 
ipjp=i j+imax+1 

i jp=i j+imaLX 

ijm=i j -imax 
im=i-l 
ip=i+l 
jm= j -1 
jp=j+l 
c 

avnx= 0 . 2 5 * ( gr adr ox { ip j ) +gr adr ox ( ip jp ) 

& +gradrox ( i jp) +gradrox ( i j ) ) 

avny = 0 . 2 5 * { gr adroy ( ip j ) +gr adr oy { ip j p ) 

& +gr adroy ( i jp) +gradroy ( i j > ) 

vx=0 . 25* (vnormx ( ipj ) +vnormx ( ipjp) 

& +vnormx ( ijp) •‘•vnormx { ij ) ) 

vy=0 .25* (vnormy ( ipj ) +vnormy (ipjp) 

& +vnormy { i jp) +vnormy (i j ) ) 

c 
c 

c.... compute the cell centered surface tension force 
c 

tens in= sigma* kappa ( i j ) /rhof 
tensx ( i j ) =avnx* tensin 

tensx ( i j ) =cvmgt (zero, tensx (ij ) ,ac(ij) . It.em6) 
c 

tensy (ij ) =avny* tensin 

tensytij ) =cvmgt ( zero, tensy (ij) ,ac(ij) . lt.emfi) 
fsv(ij ) =sqrt (tensx (ij ) **2 + tensy { i j ) * *2 ) 
c 
c 

c.... compute the cell centered shear force 
c 

f taux ( i j ) =- ( sheart) *vy/rhof 

ftaux( ij ) =cvmgt(zero, ftaux(ij ) , ac (ij ) .lt.em6) 
ftauy ( ij ) ={ sheart) *vx/rhof 

ftauy (ij ) =cvmgt (zero, ftauy (ij ) ,ac (ij ) . It .em6) 
c ftauy (ij ) =0 . OdO 

100 continue 
110 continue 
c 

i jmax=ismax (nxy , fsv, 1) 
fmax=f sv { i jmax) 
f f loor=0 . 001* f max 
do 200 j = 2 , jml 
do 200 i=2 , iml 

kappa ( i j ) =cvmgt (kappa (ij ) , zero, fsv (i j ) .gt . f floor) 
200 continue 
c 

9999 return 
end 
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but not necessarily negligible in some cases. 
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